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Abstract. We study the effects of finite temperature on the dynamics of non-planar vortices in the classical, 
two-dimensional anisotropic Heisenberg model with XY- or easy-plane symmetry. Fo this end, we analyze a 
generalized Landau-Lifshitz equation including additive white noise and Gilbert damping. Using a collective 
variable theory with no adjustable parameters we derive an equation of motion for the vortices with 
stochastic forces which are shown to represent white noise with an effective diffusion constant linearly 
dependent on temperature. We solve these stochastic equations of motion by means of a Green's function 
formalism and obtain the mean vortex trajectory and its variance. We find a non-standard time dependence 
for the variance of the components perpendicular to the driving force. We compare the analytical results 
with Langevin dynamics simulations and find a good agreement up to temperatures of the order of 25% 
of the Kosterlitz-Thouless transition temperature. Finally, we discuss the reasons why our approach is 
not appropriate for higher temperatures as well as the discreteness effects observed in the numerical 
simulations. 

PACS. 05.40.+j Fluctuation phenomena, random processes, and Brownian motion - 75.10.Hk Classical 
spin models - 75.30.-m Intrinsic properties of magnetically ordered materials 



1 Introduction 

In the past two decades, solitons and other nonlinear co- 
herent excitations have become a very generic and useful 
paradigm for intrinsically nonlinear phenomena in many 
different fields |l|0,H]. These excitations are especially im- 
portant in low dimensional systems, in terms of their re- 
lationship to key questions such as the existence of long 
range order, the mechanisms of phase transitions or the 
response to external influences . Unfortunately, for most 
problems of interest or in applications, it is not possible 
to develop an exact theory of soliton dynamics or statisti- 
cal mechanics, either because the corresponding equation 
of motion is not integrable or because perturbation terms 
added to it in order to account for relevant effects destroy 
integrability. As a consequence, much effort has been de- 
voted to develop approximate techniques allowing one to 
gain insight into soliton behavior. Among those, a very 
useful procedure is that of collective variables or coordi- 
nates H , which yields very accurate results for soliton- like 



objects with a well localized spatial structure. Collective 
coordinate techniques have, in addition, the advantages 
of their mathematical simplicity and their applicability to 
very many perturbed soliton-bearing equations, including 
most of those which are physically relevant. The validity 
of this kind of calculations has provided grounds to what 
is nowadays called the "particle-like picture" of solitons: 
In view of the fact that a global coordinate, such as their 
center of mass, is enough to describe their behavior under 
perturbations, it has been concluded that solitons can be 
treated as point-like particles in many situations. 

One important context where the above ideas are rele- 
vant is that of two-dimensional (2D) magnets and their 
collective excitations such as vortices or domain walls. 
This is a far from an academic subject: Indeed, in the last 
few years several classes of materials have been found or 
fabricated for which magnetic interactions within planes 
of their crystalline structure are much stronger than be- 
tween these planes, and therefore the magnetic properties 
are basically 2D. Materials in these classes include, for 
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instance, layered magnets (such as R^CrCU), graphite 
intercalated compounds (such as C0CI2), magnetic lipid 
layers (such as manganese stearate), and high T c super- 
conductors (see references in, e. g., ||). Many of these 
systems can be described by the classical 2D anisotropic 
Heisenberg model with XY- or easy-plane symmetry, given 

by 

H = - J E i 15 ? 1 ® + S v S v + C 1 - 5 ) S TS:} , (1) 

(m.n) 

where the subindices x,y or z stand for the spin compo- 
nents, < S < 1, and (m, n) labels nearest neighbors of 
a square lattice. Among its excitations, specially interest- 
ing ones are vortices, that are planar (i. e., with null S z 
components) if 5 > 0.297 and non-planar (i. e., with lo- 
calized S z structure) if 6 < 0.297 Such non-planar 
vortices are the specific object of our study as reported in 
the remainder of the paper; however, the ideas we will be 
discussing are general enough to be of interest in other, 
related contexts where the system behavior is governed by 
soliton-like structures. 

The first application of a collective variable technique 
to the motion of magnetic vortices and other nonlinear 
magnetic excitations was carried out by Thiele ■ For 
steady state motion, when the shape S(r, i) of the excita- 
tion in the continuum limit is rigid, he used the travelling 
wave ansatz S(r, t) = S(r — X(t)) with constant velocity 
X (the dot stands for derivative with respect to time) and 
derived the following equation of motion, 



G y x X + F = 0, 



(2) 



where F is a static force, due to either an external field 
or the interactions with other excitations. The gyrovec- 
tor Gy, in turn, is an intrinsic quantity, produced by the 
excitation itself and depending on its specific type. Gy 
is perpendicular to the XY-plane; therefore, the gyrocou- 
pling force Gy x X is formally equivalent to the Lorentz 
force. Interestingly, Thiele's equation, Eq. (|J), is first or- 
der, thus leading to non-Newtonian vortex dynamics. This 
is somewhat unusual, as in many cases solitons are found 
to behave as Newtonian point-like particles ||, obeying 
Newton's second law or its relativistic generalization. We 
return to this point below. 

The next step beyond Thiele's approach was not taken 
until very recently, when Mertens et al. jllj] developed a 
generalized collective variable theory for nonlinear coher- 
ent excitations in classical systems with arbitrary Hamil- 
tonians. Previously, Wysin et al. (l2) had tried to gen- 
eralize Thiele's equation by allowing the vortex shape to 
depend on the vortex velocity. In this way they derived 
a second order (Newtonian) equation of motion, but it 
was found that it did not agree with the simulations jyj 
U. Therefore, in fl) it was proposed that the dynamics 
of a single excitation is governed by a hierarchy of equa- 
tions of motion for the excitation center X(t). In addition, 
the Newtonian or non-Newtonian character of the equa- 
tion of motion was clarified: It was found that the type 
of the excitation determines on which levels the hierarchy 



can be truncated consistently: So-called gyrotropic (with 
|Gy| 7^ 0) excitations are governed by odd-order equa- 
tions and thus do not have Newtonian dynamics. Non- 
gyrotropic excitations are described by even-order equa- 
tions, i. e. by Newton's equation in the first approxima- 
tion. This is the situation for, e. g., domain walls. The the- 
ory in [ fj"l| was applied to non-planar (gyrotropic) vortices 
of the 2D anisotropic Heisenberg model, and it was shown 
that their dynamics is fully captured by the third-order 
equation, fifth-order corrections being negligible pT| . 

Since the zero temperature dynamics of non-planar 
vortices is completely understood, in this paper we now 
concern ourselves with the study of non-planar vortex 
dynamics in the 2D Heisenberg ferromagnet at nonzero 
temperatures. The purpose of this research is twofold: 
From the theoretical point of view, it is important to learn 
whether and when the vortex motion description in terms 
of an effective particle dynamics holds, and what are its 
main characteristics. In addition, the non-Newtonian char- 
acter of non-planar vortices could be modified by temper- 
ature, or the details of the dynamics could change as to 
eliminate the need to go beyond a first order equation. We 
note that if a collective coordinate theory at finite tem- 
perature could be worked out, it would provide a first step 
towards a statistical mechanics description of the model 
behavior in terms of a vortex gas |l4| ], as in the case of 
one-dimensional soliton bearing systems jl5| . On the other 
hand, from the experimental point of view, insofar as the 
motion of vortices has measurable consequences in inelas- 
tic neutron scattering |[(| and nuclear magnetic resonance 
experiments |l7j ] , the effects of finite temperature on vor- 
tex dynamics can have signatures in those measurements. 
The study we carry out here is then necessary if there is 
hope to compare the theoretical results to actual experi- 
ments. 

The presentation of our results proceeds as follows: 
Section | contains the study of the free and the damped 
vortex dynamics and the derivation of the corresponding 
collective coordinate theory. At this point the study is 
still deterministic, i. e., at zero temperature. Section ^dis- 
cusses how we incorporate the Langevin noise term to the 
equations of motion. Afterwards, the collective coordinate 
calculation is extended to the resulting Langevin-Landau- 
Lifshitz equation, and the mean vortex trajectory and its 
variance are computed. Section ^ contains a thorough dis- 
cussion of the comparison of the theory to the numerical 
Langevin dynamics simulation and the discussion of the 
main features of the vortex motion. Finally, Section ^| is 
devoted to the summary of our main conclusions. 



2 Zero temperature dynamics 

Our starting point is the damped Landau-Lifshitz equa- 
tion, which reads 



dS 1 



dt 



= -S Tl 



OH dS r 
x — eS x 



dt 



(3) 



where S m is the spin vector at lattice site m, H is the 
Hamiltonian, in our case that of the anisotropic Heisen- 
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berg mode l ([D ), and e is the damping parameter. Following 
Refs. [§, fllQ] and mm we have chosen Gilbert damping 
|pj| , chiefly because it is isotropic, meaning that all the 
spin components are equally damped, in contrast to the 
Landau-Lifshitz damping |2(J. As stated in the Introduc- 
tion, our approach to the problem of vortex dynamics will 
be both analytical and numerical: We first derive equa- 
tions of motion for the vortex center X(i), and afterwards 
we compare with numerical simulations for our model, i. 
e., with results from numerical integration of (^) includ- 
ing noise (see Section 0). The study of the deterministic 
(i.e., zero temperature) case we present in this section is 
an obviously necessary first step in order to be able to un- 
derstand later the problem for the full Langevin equation. 

To proceed, following [|llj we assume that the shape of 
a collective excitation depends on the velocity X and, as 
shown in , in general also on higher order derivatives of 
X(£). The corresponding generalized travelling wave An- 
satz is 



S(r, t) = S(r - X, X, X, . . . , X (n) ), 



(4) 



which yields an (n + l)-th order differential equation for 
X(t). As mentioned above, for gyrotropic excitations only 
odd-order equations are relevant, and, in the case of the 
non-planar vortices, it turned out that the third-order equ- 
ation is sufficient to describe accurately all simulations 
without damping jllj]. Therefore, in this paper we use the 
Ansatz (||) with n = 2 and apply it to the general case, i.e., 
in the presence of damping. Instead of using the Hamil- 
tonian procedure described in |jTlf , we will obtain the col- 
lective variable equations of motion in a much more direct 
way by performing the following operations with ([}]): leav- 
ing out damping for the moment, we calculate 



dS dS\ 
~dX l X ~dt) 



\dXi 

= -s' 



S x 



SH 
~5S 



, SH dS 
~5S'dX l 



dU 



= - s2 m < 5 > 



with i = 1, 2 in the case of our 2D system. H. is the Hamil- 
tonian density. According to our ansatz we insert on the 
l.h.s. 



dS 



dS 
I)X~ 



X; 



OS 



X, 



OS 



-x, 



(0) 



integrate over r and divide by S 2 . In this way we obtain 
the same third-order equation as that obtained in Ref. 
|pi[, which used Hamilton equations: 



AX + MX + GX = F 

with force F given by 



F, = - d 2 r 



on 

dXi 



(7) 



(8) 



gyrotensor G expressed as 

q. . — c-2 [ d 2 rS dS OS 

Gl '- S d dX dX 



2 ( d<f> d<P d<P dj> 

\ dX, dX, dX, dX, <' { ' 



mass tensor M with components 



Ma =S- 2 d 2 rS— x 

' dX, QX, 



d 2 r 



dip d(f) dip 



dXi dX dX, dX, 



(10) 



and third-order gyrotensor A given by 
An = S 2 d 2 rS—r- x 



dXi OX, 



= d 



dip dip d<p dip 



dXi dXi DXi dX, 



(11) 



The classical spin is constrained to have a fixed mag- 
nitude which we set to unity. Therefore, we will evaluate 
below the expressions on the right hand sides using canon- 
ical fields (p = arctan(S , y/5 , a; ) and ip = S z for the spin 
vector: 



S = yl — ip 2 cos <f> e x + yl — ip 2 sin <pe y + ipe z 



(12) 



At this time, we move on to consider the Gilbert damping 
term in (0). The same operation sequence as above yields: 



dS ( dS 
- — x S x — 

dXi V dt 



= eS' 



dS dS 
dX~~dt 



eS 2 



dS dS ■ 
- X 



dS dS •• 

— Xj 



dXidXj J dx.dx 



dS dS 

dXdX. 



■A, 



J(13) 



An integration over r gives three terms which can be com- 
bined with the three terms on the l.h.s. of (Q), i. e., the 
damping appears in every order 

(A + a)X + (M + m) X + (G + g) X = 

= AX + MX + GX = F. (14) 

The components of the damping contribution to the ten- 
sors are 



9ij = e Jd 2 



dS dS 

dXidXj 

o. ld 2 r \ (l-^ 2 ) 6 " 



1 



dip dtp 



dX t dXj 1 - ip 2 dX, dX t 



(15) 



= e d 



dS dS 

dXidXj 

< ,d 2 r{(l-ip 2 ) d(b 1 



dip dip 



dX z dX 



1 - ip 2 dX, dX t 



(16) 



dfr- 



dS dS 

dXdX 



dip dip 



dX, dX 1 - ip 2 dXi dX t 



(17) 
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We note that the first-order part of (14) was already de- 
rived by Thiele B. 

Now, we address the problem of the explicit calculation 
of all the tensor components. This is possible only if the 
dynamic structure of the collective excitation is known. 
The Hamiltonian density derived from ([I]) reads 



Mi. 



ML 



4o a r 



A 



G 
16S 



Ca, 



(28) 



(29) 



n 



JS 2 



</> 2 )(V0) 2 +<5[4i/> 2 - (VVO 5 



2, L 
eirq m — 

a c 



a 



Si 



(30) 



wry (is) 



G 

m lj =me lj , m = e— (i? - a 2 ) + C m , (31) 



In ||ll| the Hamilton equations were considered for a non- 
planar vortex in the center of a circular system with free 
boundary conditions. The vortex structure is complicated 
in an inner region < r < a c « 3r v , where 



and 



dij = aSij, a = e 



L In L — a c In a c ) 



1 



(19) 



characterizes the vortex core M. 5 is the anisotropy pa- 
rameter in ([jj). Recalling that non-planar vortices are sta- 
ble for < S < 0.297 for a square lattice, we will use 
5 = 0.1 for our simulations. We note that the inner re- 
gion contributes very little to the integrals in (10), (11) 
and (15)-(17); except for (9), the dominant contributions 
stem from the outer region a c < r < L, if we choose a 
large system radius L. Here the vortex has the following 
dynamic structure, which is known to be a very accurate 
description from simulations |ll| : 



/>0 + 01 + ^ = ^0+^1+^2 (20) 



with 



9o = gtan — . 

Xi 



!>i = p{xiXi + x 2 X 2 ), 



q v 

^2 = ^7 In — (x 2 Xi - xiX 2 ), 
8o eL 



tpo ~P\j — exp(-r/r v ), 



Vi 



4<5r 2 



(x 2 Xi - X1X2), 



and 



^2 = tt( x i^i + x 2^ 2 ). 
4o 



(21) 



(22) 



(23) 



(24) 



(25) 



(26) 



Here q = ±1 is the vorticity and p = ±1 is the polar- 
ization, which determines to which side the out-of-plane 
structure of the vortex points. Straightforward integra- 
tions then yield the expressions of the tensor components: 



Gi 



Gdj, G = 2irpq, 



(27) 



i (L 2 -a 2 ) \+C a , (32) 



where 6^ is the 2D unit matrix, is the antisymmetric 
tensor, and the different constants C are the contribu- 
tions from the inner region of the vortex. We see that in 
every odd-order of (14) a symmetric damping matrix is 
combined with an antisymmetric normal (non-damping) 
matrix, and vice versa for the even orders. Moreover the 
size dependence of the n-th order damping components is 
the same as that of the (n + l)-th order normal compo- 
nents. The first-order damping elements (30) were already 
evaluated in Q and p0{ . 

For the solution of the equation of motion (14) we pro- 
ceed as in Ref. [[ll]: We consider small displacements x 
from a mean trajectory X°, on which the vortex is driven 
by F 



X(t)=X°(t)+x(t). 



(33) 



We will denote the components of x by x± and x 2l with 
the caveat that they should not be confused with the com- 
ponents of r in Eqs. ( pl| ) through (26). In view of our 
simulations, we consider the situation where the force is 
always pointing in the Xi-direction and expand to first- 
order around -Xi(O) = Rq (this is justified because in our 
simulations Fq, and even more Fq, is very small), i. e., 



F = F + Ft>x l . 



(34) 



For X®(t) we obtain two coupled linear third-order equa- 



tions. Taking initial conditions X®(0) 
the solutions are 



F Q 



X<t = R + -£[e X p(t/T)-l] 



GF 



i[exp(t/r)-l], 



Ro, X°(0) = 

(35) 
(36) 



where r is determined by a cubic equation. The mean 
trajectory is a straight line X 2 = G/g(X® — Rq), which 
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slightly deviates from the X2-axis. The angle g/G is small 
because g ~ e, where we choose small damping constants 
e in the simulations. As r is of the order of G 2 /{gF^) it is 
very large, in fact much larger than our integration times. 
Therefore one can expand ( |35| ) and one get a constant 
velocity on the mean trajectory: X® — gFo/G 2 , X° = 
F /G. 

The motion around the mean trajectory is obtained 
by solving the two coupled linear third-order equations 
for the displacements x(f ) using the Ansatz 



— x\ exp[— ((3 — iuj)t] 



(37) 



We find 



0- 



±iM- 



2{A±ia) 



± 



± 



y/{±iM + m) 2 - A{A ± ia)(G ± ig) 
2{A±ia) : 



(38) 



with amplitude ratios k = x\jx\ = ±1 and phase dif- 
ferences ±7r/2, where we have set Fq = for simplicity. 
With Fq 7^ Eq. (38) becomes even more complicated and 
| k | 7^ 1. The separation of real and imaginary parts leads 
to cumbersome formulas. Therefore we compute the fre- 
quencies Wi ( 2 and the relaxation constants /3i,2 as a func- 
tion of the parameters e and L; we choose q = p = 1 for the 
charges and <5 = 0. 1 for the anisotropy. The a c -dependent 
parts in (p8|)-(|3l|) can be combined with the constants Cm 
etc.; the combined constants can be neglected for large 
systems. As for the frequencies, wi_2 turn out to be very 
close to each other; hence, the important parameters will 
instead be their mean and difference. Examples of their 
numerical values for a system of radius L = 24 are 



= >/(Jl£J2 ~ 0.05, AuJ = LU2 — Wi » 0.01, 



(39) 



for a wide range of damping values (up to e = 0.05), 
whereas for fixed e the frequencies decrease with 1/L up 
to rather large systems (L = 5 000) . Plots of lo c and Aw as 
a function of e and L can be found in pl[ (note, however 
that the caption under Fig. 1 of must read 48 instead 
of 24). 

In the simulations the purpose of the damping is to 
dissipate the energy which is supplied to the system by the 
noise. Therefore we must know the range of e (for a given 
system size) in which the frequencies are not influenced 
by the damping. As shown in [ Ell , this range is defined by 
the condition 



eL < 5. 



(40) 



The relaxation constants are nearly equal and the 
mean value is [3 C = e/8, and for the above range w c and 
Auj are related to the parameters G, M, and A in a very 
simple way O] 



Au 



M 
T 



hxL 




(41) 



Fig. 1. Sketch of the vortex motion as governed by the Landau- 
Lifshitz equation with Gilbert damping. The plot is approxi- 
mate and does not correspond to an actual simulation. 



Finally we briefly discuss the shape of the trajectories. 
We first consider the motion in a frame which is mov- 
ing along the mean trajectory X°(i): The general solu- 
tions for the displacements Xi(t) are linear superpositions 
of ( |37| ) with uji.2- Both xi(t) exhibit a very pronounced 
beat because wx,2 are nearly equal. The orbits X\{x2) are 
Lissajous curves, which can look very intricate for certain 
parameter ranges. We go into the laboratory frame by 
adding X°(t). Without the splitting of ui\,2 we would get 
a cycloid. Due to the splitting we finally get a superpo- 
sition of two cycloids, which are damped because of 0\ t i- 
A cartoon of the vortex motion on a circular system is 
sketched in Fig. |l|. Here the vortex is driven by an im- 
age force which points in radial direction (see Section |]). 
Without damping, the mean trajectory X°(t) would be a 
circle, due to the gyrocoupling force. With damping, the 
circle converts to an outward spiral. However, the radial 
motion is very exaggerated in the sketch, and the same is 
true of the damping of the cycloidal oscillations around 
the mean trajectory. The amplitude of these oscillations 
in fact remains of the order of a lattice constant for a long 
time. 



3 Finite temperature dynamics 

3.1 Derivation of the vortex equation of motion 

In order to study the finite temperature dynamics of vor- 
tices, we must introduce thermal noise in the Landau- 
Lifshitz equation with damping, Eq. (||). However, we can- 
not simply add independent noise terms to these three 
equations, because if we do so we do not arrive to some- 
thing in the form of Langevin equations (all components of 
dS/dt appear in each equation due to the cross-product). 
Therefore we must first take all dS/dt-terms to the l.h.s., 
casting it explicitly into a first order equation, and only 
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then introduce the three white-noise terms rf a (r, t) , yield- 
ing: 



dS 



1 



1 + e 2 S 2 



S x 



5_H_ 
~5S 



eS x 



S x 



SH_ 
~5S 



V > (42) 



with 



W a (r,t))=Q t (43) 
(rf a (r,tW p (T ',t')) = 2ek B T6(r'- r)5(t' - t)S af} , (44) 

where a, p — 1, 2, 3. Now we take rj' to the l.h.s. and 
undo the above procedure, i. e., we write (42) in the same 
form as (||), but with S — rj' instead of S. We thus arrive 
at 

dS „ 5H _ dS 

d^ = - SX ^- eSX d^ +T? (45) 

with 



where summation over repeated indices is implicitly un- 
derstood. The mean is easily shown to be zero, whereas 
for the correlation functions j22| , from Eq. (Q) we obtain 



(*f(*)W)) = 

= 2ek B T5(t 







dVi Q) (r)/, 



(a) 



(r) . (51) 



Instead of the S a we introduce the fields 4> and ip in 
(O) and thereby fulfill the constraint |S| = 1. After some 
algebra we obtain 



Var(Ff) = 2efc B T J d 2 r j(l - ^ 



2^ (9± 
\dXi 



1 



■ (52) 



1- V 2 \dX, 



We note that in this equation the leading contribution 
comes from the static vortex structure, as given by Eqs. 
©and©: 



rj - e(S x n') 



(46) 



If we now compute the variances of r/, we find that the 
width of the distribution for the component parallel to the 
spin vector is <Tq = \f2ek~oT, w hile the w idths for the per- 
pendicular components are coVl + e 2 «S 2 . In the Langevin 
dynamics simulation we will apply the constraint |S| = 1, 
which means here that only the perpendicular components 
are relevant. Thus we can replace r] in (pl5| ) by T]' if we cor- 
rect the widths by a factor of ylT? . Taking into account 
that in our simulations we will be using values of e of the 
order of 1CP 3 (see Sec. Q), we will neglect the correction 
factor in the following. 

As in the previous section, we calculate 



S(^'l| = (S^]'l 



dS 



(47) 



integrate over r and combine this with the results of the 
previous section; we have thus found the collective coor- 
dinate equation in the presence of noise, namely 



AX + MX+GX = F + F st , 
where the stochastic force is given by 



(48) 



(49) 



To achieve a complete understanding of the vortex dynam- 
ics as described by Eq. (ft8l), we need to know the mean 
(Ff 1 ) and the variance Var(Fj St ). We define 



(a) 



1 dS 7 



(50) 



Var(i^ st ) = 2irek B T / dr r 



I -Mr? 



1- Mr) 2 



(53) 



As tpo decays exponentially, the second integral is inde- 
pendent of L, while the first one grows logarithmically. 
This suggests that, in order to approximately calculate 
Var(i 7 ! i st ), we can divide the integral in an outer part from 
a c < r < L and a core part C(a c ). By doing so we can 
write 



Var(Ff) = 2ek B T ■ tt 



In — + C(a c ) 



(54) 



which implies that the stochastic forces can be represented 
as white noise on the level of the collective coordinates 
with the properties (F- St ) = and 



(Ft(t)Ff(t')) = D v S ij 5(t-t% 



(55) 



where the effective vortex diffusion constant Dy is deter- 
mined by the r.h.s. of j54|). We recall that the diffusion 
constant D on the microscopic level, i. e., the one we will 
use in the simulations, is D = 2ek B T. 

The core contribution C{a c ) in the integral (53) cannot 
be calculated accurately by using ipo(r) from the contin- 
uum limit |0]. Therefore, we have computed the full inte- 
gral (53) using for ipa an ad-hoc function which was fitted 
to the static vortex structure as obtained from the simula- 
tions at zero temperature |^l|,^4). The results for L = 24 
and 8 = 0.03,0.10, and 0.30 are D v / D = 10.02,12.08, 
and 14.18, respectively. 
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3.2 Solution of the equation of motion 

We now turn to the solution of Eq. (Est). As the force F, 
which drives the vortex, can be expanded up to first order 
around the mean trajectory as discussed in Sec. ^, we find 
the linear equation 



with 



Mx + Gx fx = F s 



F>0 




(56) 



(57) 



Thus we can use the Green's function formalism to obtain 
a formal solution. This proceeds in two steps. First, the 
Green's function matrix is obtained from the solution to 
the equations 



Agj + Mg, + Gg 4 - f gi = S(t)Ii, 



(58) 



where gi, with i — 1,2, are the two columns of the Green's 
function matrix, Ij are the corresponding columns of the 
identity matrix, and suitable conditions have to be im- 
posed on both eqs. (|5q). Once the Green's matrix G has 
been calculated, the second step is to solve the stochastic 
problem (|56"|) . Its solution is then exactly given by 



x(t) =x ft (t) 



dsG(t-s) F st (s), 



(59) 



where x^(i) stands for the solution of the homogeneous 
version of Eq. (|56|). We note that G should not be confused 
with the gyrocoupling tensor in Eqs. (0), (9), and (p7|). 

Let us now discuss the first part of the calculation, i.e., 
the computation of the Green's matrix G. Eqs. ( |58| ) above 
need to be suplemented with the following conditions: G, 
G, and G vanish for t < 0, and 



G(0+) = G(0+) = 








G(0+) = A" 1 = 



1 



-A 

i 2 + A 2 \ A a 



(60) 
(61) 



In order to find the columns g,; of the Green's matrix, we 
take the Ansatz 



do not present it here insofar as the derivation is straight- 
forward. 

Once the c k and hence G are known, we can move to 
the second part of the procedure, namely to find the tra- 
jectory and to evaluate its relevant moments. It is evident 
from Eq. (|59|) that the mean trajectory will be exactly 
the same as that of the deterministic case, because the 
average of the integral of F s * vanishes. We will therefore 
concentrate on the variances, 



At) = {xiXj) - (x l )(x j ). 



(63) 



Using the expression ( p9[) it can immediately be seen that 



;(*) 



2 

E 

k=l 



dt'D v G ik {t-t')G jk (t-t'), 



(64) 



where Gij stand for the elements of the Green's matrix. 
Once again, the calculation is simple but tedious, due to 
the many terms involved by the product of the Green's 
matrix elements. Aside from this, the expression is easily 
obtained as the integrals involve only exponentials. As an 
example, we present a summary of the calculation of a\ x , 
which is the simplest element of the variance matrix. Nev- 
ertheless, in order to facilitate the presentation and the 
subsequent discussion we have made the following simpli- 
fications: i) u)\ = u>2 = u c and 0\ = 02 = Pc because 
the splittings are very small [see Eq. (^) and above]; ii) 
lu 2 + 2 ~ u 2 because C = e/8 and e = 0.002 in the sim- 
ulations, implying C is two orders of magnitude smaller 
than ui c as given by Eq. (|39"|); iii) A 2 + a 2 ~ A 2 , because 
a/A = O(e), see Eqs. ( |29| ) and (32), and iv) exponential 
terms involving tjr are expanded to first order, because r 
is much larger than our integration times; see below Eq. 
36[). Within these approximations, it can be shown that 



°ii(t) 



1 

4ft 



(1 



-20 C U 



-Wet 



4w c 



sin 2cj r t 



(65) 



For small times, t <§C 1/0 C > we are left with an expression 
which implies linear behavior plus oscillations, given by 



Si(t) 



6 

£< 

k=l 



exp(A fe i) 6{t), 



(62) 



where 9(t) is the Heaviside function, at and b k form the 
eigenvectors belonging to the eigenvalues of the homo- 
geneous problem, i.e., Eq. (bq) with its r.h.s. set to zero, 

(i) — 

and c k arc the unknown amplitudes in the linear combi- 
nation. The eigenvalues are already known from the previ- 
ous section: one is zero, one is 1/r [see below Eq. (|36|)], and 
the other four are give n by Eq. (38). All that remains is to 
insert the Ansatz (J62J) in Eq. (|5q) and find the values for 

from the corresponding system of algebraic equations. 
Their expression is rather cumbersome and therefore we 



(*) = 



Dyt 

Ah4 



3 sin io r t 

77 - 2 ~ 

2 u> r t 



sin 2w c t 

40J r t 



(66) 



We note that this function increases monotonously and 
that it has no extrema but only inflection points. For large 
times, t 3> l/0c, only the first term of Eq. (65) remains, 
and the variance becomes a straight line, 



°Ut) 



Dyt 



Dy 

G 2 



t. 



(67) 



where Eq. ( fi"l| ) has been taken into account. Interestingly, 
this result is identical to the one obtained by omitting the 
second- and third-order terms in the vortex equation of 
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motion (|48|). We thus conclude that these terms have two 
effects: First, they produce the oscillatory parts in Eq. (65) 
(note that they are naturally connected to the cycloidal 
vortex trajectories), and second, for small times the slope 
of a\ x in Eq. (|66|), averaged over the oscillations, is larger 
by a factor of 3/2 compared to Eq. (p7|). 

We do not present here the other elements of the vari- 
ance matrix because they contain even more terms than 
Eq. (65). Instead, we simply record the expression for their 
long time behavior, which is 



'12 



1 22 



D V F .2 

G 2 2G ' 



Dv 

Q2 



t + 



§) r 



(68) 
(69) 



The quadratic and cubic terms in t appear in addition to 
the standard random walk result which is proportional to 
t. These additional terms arise because we have allowed 
that the driving force F depends on the vortex position, 
see Eq. (HJ). We have considered a force in the X\ direc- 
tion which drives the vortex in the X2 direction, due to 
the gyrocoupling force Gy x X in Eq. (||) or GX in Eq. 
(Q), respectively. Therefore, only the 2-components of a 2 
are affected, a 2 2 acquiring a factor (Fq/G) t, 022 acquiring 
it twice. 



4 Langevin dynamics simulations 
4.1 Numerical procedure 

We begin with one vortex with its center located at a dis- 
tance Rq from the middle of a circularly shaped square 
lattice with a radius of L lattice constants. We use free 
boundary conditions to get an image antivortex which 
leads to a radial force on our vortex (see Jll|,^| and be- 
low). The initial spin configuration stems from an itera- 
tive program which produces a discrete vortex structure 
on the lattice [^4|. In this way we avoid the radiation of 
spin waves which would appear during the early time units 
if we use a continuum approximation for the vortex shape. 
The parameter ranges must be chosen very carefully for 
the following reasons: i) We want that the vortex moves 
smoothly over the Peierls-Nabarro potential of the lattice; 
hence, the diameter 2r v of the out-of-plane structure must 
be considerably larger than the lattice constant. Setting 

5 = 0.1 we find 2r v ~ 3 from Eq. (|l9|); ii) we choose a 
system radius L = 24 which provides enough space: the 
vortex moves outwards roughly on a spiral, but even for 
very long integration times the out-of-plane vortex struc- 
ture should not contact the boundary, and iii) for the same 
reason the initial distance Ro from the middle of the cir- 
cle should not be too large. On the other hand, Ro should 
not be too small; otherwise the driving force F would not 
be strong enough to overcome the pinning forces of the 
lattice. Letting Ro — 10 both conditions can be simul- 
taneously fulfilled, if the damping e is small enough (the 
larger e, the sooner the vortex reaches the boundary). Note 



however that a small e means long saturation times (see 
below) . 

For the time integration of the Landau-Lifshitz equa- 
tion we use the discrete version of (42) where dS/dt has 
already been isolated on the l.h.s. . In contrast to our an- 
alytical calculations we work here with the cartesian com- 
ponents S a . Therefore we explicitly take into account the 
constraint S 2 = 1 by adding S times a Lagrange param- 
eter A to (42), see Ref. J2j|. We form the time derivative 
of the constraint, eliminate A and get 



with 



U 



e 2 S 2 



dt 



-S x 



U 



SH_ 

Is 



su, 

s 2 ' 



eS x 



(70) 



S x 



6H_ 

Is 



(71) 



where the site index has been omitted. We note that ( |7(i| ) 
is the same as orthogonalizing S and U by the Gram- 
Schmidt method. For the time integrations we use the 
same code as in |]lTf . In addition, the position of the vortex 
center, in particular the position within a lattice cell, is 
determined by a procedure also discussed in [ll|. 

To find a proper damping constant we checked the time 
dependence of the system energy using different damping 
constants for L = 24 and T = 0.02. The energy at t = is 
the same as for T — and e = because the noise will be 
introduced with the first time step of the simulation. The 
energy then rises and saturates on a value independent of 
e, but for e > 8 x 10" 3 the energy decreases slowly after 
saturation. The saturation time gets longer with lower e, 
for e > 2 x 10~ 3 we achieve acceptable saturation times 
< 300 (in units of h/(JS)). We have always made a pre- 
run of length to > 300 prior to beginning the evaluation 
of the simulation data. 

The difference between the energy without tempera- 
ture and the saturation energy with temperature must be 
the thermal energy. We computed the mean thermal en- 
ergy per spin at several temperatures and it agreed with 
J/2 x fc B T up to T = 0.9, / being the number of degrees 
of freedom per spin. For higher temperatures we find too 
low values for the energy. We believe that the numerical 
procedure would have to be improved if we were interested 
in this regime. 



4.2 Vortex trajectory 

We studied the trajectory of the vortex center at different 
temperatures keeping L = 24 and e = 2x 10~ 3 fixed. We 
can distinguish three temperature regimes in which the 
trajectories differ qualitatively. 

For < T < T3 as 0.05 we observe two frequencies 
in the oscillations around the mean trajectory which can 
be identified with the cycloidal frequencies lu±.2 in (38). 
The intensities of the Fourier peaks at uj\£ decrease with 
temperature and vanish at T3 in the background, but lo\^ 
are constant in the whole regime. This means that here 
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Fig. 2. Average trajectory of a vortex for temperature T = 
0.03, damping e = 0.002, system radius L = 24 and an ensem- 
ble of 100 realizations. Upper panel: Radial coordinate of the 
vortex center vs time. Middle panel: Azimuthal displacement 
ip = <f)(t) — ujot, where loo is the angular velocity on the mean 
trajectory of radius Ro- Lower panel: Fourier spectrum of R(t) 
in the upper panel. The spectrum of ip(t) is very similar. 



the third order equation of motion ( |48| ) with temperature 
independent parameters can describe the vortex dynam- 
ics. For one temperature of this regime we plot in Fig. 
U the average radial coordinate R(t) and the azimuthal 
displacement <p(t) = </>(i) — u>ot. We want to stress that 
the plots present averaged results: For the computations 
of the vortex trajectories and variances we have averaged 
over 100 different runs starting from the same initial con- 
figuration, which is defined as the final configuration after 
a pre-run of length to = 1250. In the expression for cp(t) 
above, loq = Fq/{GRq) is the frequency of the rotation on 



the mean trajectory which is essentially a circle where the 
radius Rq grows very slowly with rate gF /G 2 due to the 
damping. On the mean trajectory the vortex is driven by a 
radial force Fq due to the image vortex at i?W = L 2 /Rq, 
which has opposite vorticity but the same polarization 
[ p6[ . As the average motion is very slow (uj rs 2.5 x 10~ 3 ) 
we can actually work in a cartesian system and use the 
results (|33|)-(38). Here the ATi-axis points in the radial di- 
rection, and the A^-axis in the azimuthal direction |j27f . 
The lowest panel of Fig. ^| shows the Fourier spectrum of 
R(t). In addition to U12 one also observes the difference 
Aui = LU2 — uj± . This can be explained by working in polar 
coordinates, which is not discussed here because the for- 
mulas become much too cumbersome. The peaks at higher 
frequencies are second harmonics of LO\p- 

For T3 < T < T\ as 0.3 we do not observe the above 
mentioned two frequencies any longer. In this regime, we 
found that some runs had to be excluded from the aver- 
age because the vortex suddenly changed its direction of 
motion. This occurs because, opposite to the case of the 
vorticity g, the polarization p of the vortex is not a con- 
stant of motion for a discrete system: The out-of-planc 
vortex structure can flip to the other side of the lattice 
plane due to the stochastic forces acting on the spins. 
Then G = 2irqp in Eq. (g^) changes sign and thus the 
direction of the gyrovector in Eq. (^|) is reversed, which 
implies that the direction of the vortex motion is reversed 
as well. This noise-induced switching between the two vor- 
tex polarizations is a very novel effect in itself, and hence 
we are developing a theory for the switching rate psfl . In 
this respect, it can be mentioned that switching can also 
be induced by an ac magnetic field in the easy plane. As 
the symmetry is broken here, such a switching occurs only 
for one sense of rotation, and there is no transition back 
to the original state ]2jJ. 

Finally, for T > T\ , a single- vortex theory as presented 
here is no longer adequate because at these temperatures 
the probability for the spontaneous appearance of vortex- 
antivortex pairs becomes too large. These pairs can break 
up above the Kosterlitz-Thouless transition temperature 
Tkt — 0.7 in our units. Between T\ and Tkt, these pairs 
interact with the initial vortex although they are not sep- 
arated, thus introducing new forces and effects which the 
present theory does not take into account. Moreover, very 
recent Monte Carlo simulations 



have revealed 



31 



that for higher temperatures the vortex motion is strongly 
influenced by creation and annihilation processes: Typi- 
cally, an unbound vortex travels only one or a few lattice 
spacings until it annihilates with the antivortex of a pair 
which meanwhile appeared spontaneously in the neighbor- 
hood. Then, the vortex of this pair continues the travel 
instead of the original vortex. 



4.3 Variances of the vortex trajectories 

As the vortex positions in the simulations are evaluated 
in polar coordinates, we obtain a variance matrix with 
a% A and cr? . . Their time evolution is plot- 



elements Crfy 



'Rcb 



ted in Fig. H for T = 0.03, which is close to the upper 
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Fig. 3. Variances of the vortex trajectory; parameters are 
the same of Fig. |^. From top to bottom, shown are the 
variance of the radial coordinate, <j\ r = {R 2 } — (R) 2 
the off-diagonal elements of the variance matrix, and a, 
{(Ro<p) 2 ) - {{Ro<p)) 2 - In all three cases the lower line is the 
theoretical prediction without adjustment of parameters. 



edge of the low temperature regime defined in the previ- 
ous subsection. The solid lines are the theoretical results 
from Section 0, without the simplifications i)— iv) discussed 
above Eq. (65), which were only made there to facilitate 
the discussion. As the theory has been worked out in carte- 
sian coordinates, the following factors appear when go- 
ing over to polar coordinates: No factor in o 2 RR , a factor 
K = 1 — Fq/(FqRq) in cr^, , and a factor k 2 in front of the 
terms cubic in time in a"i + . 

Figure || shows that, for not too long times, the agree- 
ment between theory and simulation is astonishingly good; 
it is important to stress that no parameters were adjusted 



at all. Moreover, we worked in the continuum limit, while 
the simulations were performed on a discrete system. For 
very long times, t > 2 000, the agreement becomes poorer. 
This is partially due to one simplification of the theory, 
namely that we have used a constant Rq although during 
the simulation Rq slowly increases by several lattice con- 
stants as the trajectory is roughly a spiral (see the cartoon 
in Fig. |l|). The force term Fq — F'(Rq) increases as well, 
because the force increases when the distance to the image 
vortex becomes smaller. As Fq appears in crf^ and (Fq) 2 
arises in front of the cubic term in a 2 ^ in Eqs. (|6^) and 
(|69|), including this effect would lead to an improvement 
of the agreement between theory and simulations. 

Aside from those discussed above, there is another pos- 
sible reason for the discrepancy between theory and sim- 
ulation whose consideration, unfortunately, would lead to 
very involved calculations: The integral (52) for Dy [as 
well as the integrals in Sec. || except (9)] have been evalu- 
ated by placing the vortex into the middle of the circular 
system. However, in the simulations the distance from the 
lattice center is Rq, which moreover increases slowly. We 
have estimated the above integrals by expanding in Rq/L, 
which shows that the first order terms vanish. Neverthe- 
less, the second order terms yield corrections which are al- 
ready of the order of 20% for Rq — 10, becoming larger as 
Ro increases. Even more, the variance (|54|) of the stochas- 
tic forces is actually a diagonal tensor, see Eq. (51) and 

[ pT| | . Therefore, we get a radial diffusion constant Dy 

which differs from the azimuthal constant Dy ^ when the 
vortex is not at the center. This splitting is also of order 
(Rq/L) 2 . As Dy ' appears, e. g., in front of the cubic term 

in cr^, whereas the linear term contains Dy , it is quite 
possible that the agreement with the simulations could 
be improved by taking into account the splitting of the 
diagonal elements of the diffusion tensor. 

We numerically integrated up to times t = 4 000 (let us 
point out in this regard that this takes three weeks CPU 
time on a CRAY-YMP/EL for averages over 100 runs) be- 
cause this is the characteristic time given by 1 / ' (3 C — 8/e for 
the damping in the trajectories. We should see then that 
the slope of the time-averaged function o- 2 RR gradually de- 
creases, eventually by a factor one third for t ^> 1/(3 C [cf. 
the discussion of Eq. (65)]. We checked this for the theo- 
retical results in Fig. |3J, and found that in the simulation 
data this effect can be observed only qualitatively. For <t|^ 
and the effect is hidden by the quadratic and cubic 
terms. 

We would like to stress that the strong fluctuations in 
Fig. H (which seem to be smaller in the two lower plots 
because of the different scales) arise not only due to the 
noise but also from discreteness effects. This is demon- 
strated very clearly by the simulation in Fig. || for a very 
low temperature (T = 0.003) using 1 500 realizations. We 
have identified the sharp spikes as discreteness effects by 
comparing with the times when the vortex center moves 
over the ridges along the lattice lines (these times are indi- 
cated as dashed vertical lines) . The vortex energy is high- 
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Fig. 4. Variance of the radial vortex coordinate averaged over 
1 500 realizations, for T = 0.003, e = 0.002 and L = 24. The 
dashed vertical lines indicate the times at which the vortex 
center moves over ridges of the Peierls-Nabarro periodic po- 
tential. 
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Fig. 5. Vortex diffusion constant Dv as a function of tem- 
perature, for e = 0.002 and L — 24. Solid line: Theoretical 
results from Eq. (53); dashed line: Adjusted Dv from fitting 
the theoretical curves for a 2 (t) to the simulation data. 

est when the center is at a lattice point, and lowest in the 
middle of a cell. 

Last, but not least, we discuss the temperature depen- 
dence of the vortex diffusion constant Dy ■ A linear depen- 
dence is predicted by Eqs. (53) and (|54|). For comparison 
with the simulations we have fitted the theoretical curves 
to the observed variances by adjusting Dy, which appears 
as a factor in front of all the components of a 1 . This was 
done for two temperature decades. Fig. || shows a nearly 
linear dependence, and therefore the only difference be- 
tween Dy from the simulations and the theoretical Dy is 
a constant factor of about 1.8 for the whole temperature 
regime. 



5 Conclusions 

In this paper, we have reported our analytical and numer- 
ical work regarding the effects of temperature on the dy- 
namics of non-planar vortices in 2D, classical, anisotropic 
Heisenberg ferromagnets. As a preliminary result, we have 



described the zero temperature dynamics of vortices in the 
presence of Gilbert damping. We found that damping con- 
tributes to all the terms of the third order equation of mo- 
tion for the vortex position, but its contribution is always 
an order smaller in the system size than the correspond- 
ing free propagation part. We have solved the equations 
of motion and qualitatively discussed the motion of the 
vortex, which consists of a mean straight trajectory plus 
(damped) additional oscillations. By means of the same 
analytical approach, we have been able to derive a third 
order stochastic equation of motion for the vortex center 
when thermal noise is added to the system. The equation 
shows that the effective stochastic force acting on the vor- 
tex is also a Gaussian white noise, whose variance depends 
linearly on the temperature. We have exactly solved the 
stochastic equation of motion and obtained analytical ex- 
pressions for the mean vortex trajectory and its variance. 
The variance along the coordinate perpendicular to the 
direction of motion of the vortex is diffusive, i.e., it in- 
creases linearly with time; however, other components of 
the variance matrix (the parallel-perpendicular and the 
parallel-parallel terms) turn out to include nonlinear con- 
tributions coming from the fact that the vortex motion is 
perpendicular to the driving force, due to a Lorentz-like 
gyrocoupling force. 

The above summarized analytical results, obtained in 
the continuum limit of the Landau-Lifshitz equations gov- 
erning the model dynamics, have been compared to Lan- 
gevin dynamic simulations of the discrete 2D Heisenberg 
model. The numerical results allow us to establish three 
different temperature regimes for the vortex propagation: 
a low temperature one, where the vortex motion follows 
essentially the third order equation of motion with pa- 
rameters independent of temperature; a middle tempera- 
ture one, at which traces of the oscillations arising from 
the third order equation are lost, and a high temperature 
regime, which is not describable by a one- vortex approach 
because too many vortex-antivortex pairs arise in the sys- 
tem. Our analytical results are seen to be a good descrip- 
tion of the vortex dynamics up to temperatures of the 
order of 10% of the Kosterlitz-Thouless transition tem- 
perature. Remarkably, the analytical predictions, which 
include no adjustable parameters, agree qualitatively well 
with the numerical simulations, and even quantitatively 
at early times. The agreement becomes worse for longer 
times due to the approximations involved in our theory: 
The calculations were made for a constant radius of the 
trajectory and a constant force F gradient, aside from sim- 
plifications necessary to calculate the integrals which give 
the parameters for the equation of motion. In addition, we 
have been able to clearly identify the influence of discrete- 
ness in the numerical results, which cannot be captured by 
our continuum theory. Finally, we have also verified that 
the vortex diffusion constant depends linearly on temper- 
ature as predicted, although the quantitative comparison 
is not correct by a factor two. We thus conclude that the 
collective coordinate theory we have derived for vortex dy- 
namics is a good description of the phenomena observed 
numerically at low and intermediate temperatures. The 
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discrepancies between theory and simulations have been 
understood in terms of the unavoidable approximations 
involved in the calculations. Finally, we note that for vor- 
tices quantum effects are possibly more important than 
for kinks in one-dimensional spin models, where at least a 
part of these effects can be taken into account by a renor- 
malization of the kink parameters. For 2D spin models, it 
is not clear how a quantum vortex should be defined. In 
any case, a finite lifetime and other novel features seem to 
appear p3[ . 

We thank Esteban Moro and Grant Lythe for discussions. Tra- 
vel between Bayreuth and Madrid is supported by "Acciones 
Integradas Hispano- Alemanas" , a joint program of DAAD (Az. 
314- AI) and DGES. Travel between Europe and Los Alamos is 
supported by NATO grant CRG 971090. Work at Madrid and 
Leganes is supported by CICyT (Spain) grant MAT95-0325 
and by DGES (Spain) grant PB96-0119. Work at Los Alamos 
is supported by the United States Department of Energy. 



References 

1. Nonlinearity in Condensed Matter, edited by A. R. Bishop, 
R. Ecke, and J. Gubernatis (Springer, Berlin, 1993). 

2. Nonlinear Coherent Structures in Physics and Biology, 
edited by K. H. Spatschek and F. G. Mertens (Plenum, 
New York, 1994). 

3. Fluctuation Phenomena: Disorder and Nonlinearity, 
edited by A. R. Bishop, S. Jimenez, and L. Vazquez (World 
Scientific, Singapore, 1995). 

4. See, e.g., A. R. Bishop, in Ref. @ 

5. A. Sanchez and A. R. Bishop, SIAM Review in press 
(1998). 

6. A. R. Volkel, F. G. Mertens, A. R. Bishop and G. M. Wysin, 
Phys. Rev. B43, 5992 (1991). 

7. M.E.Gouvea, G. M. Wysin, A. R. Bishop and F. G. Mer- 
tens, Phys. Rev. B39, 11840 (1989). 

8. G. M. Wysin, Phys. Lett. A, in press. 

9. A.A.Thiele, Phys. Rev. Lett. 30, 230 (1973). 

10. A.A.Thiele, J. Appl. Phys. 45, 377 (1974). 

11. F. G. Mertens, H.-J. Schnitzer and A. R. Bishop, Phys. 
Rev. B 56, 2510 (1997). 

12. G. M. Wysin, F. G. Mertens, A. R. Volkel and A. R. Bishop, 
in §. 

13. A. R. Volkel, G. M. Wysin, F. G. Mertens, A. R. Bishop, 
and H. J. Schnitzer, Phys. Rev. B 50, 12 711 (1994). 

14. F. G. Mertens, A. R. Bishop, G. M. Wysin, and C. Kawa- 
bata, Phys. Rev. B 39, 591 (1989). 

15. J. F. Currie, J. A. Krumhansl, A. R. Bishop and S. E. 
Tmllinger, Phys. Rev. B 22, 477 (1980). 

16. K. Hirakawa, H. Yoshizawa, J. D. Axe, and G. Shirane, 
Suppl. J. Phys. Soc. Jpn. 52, 19 (1983); L. P. Regnault, J. 
P. Boucher, J. Rossat-Mignod, J. Bouillot, R. Pynn, J. Y. 
Henry, and J. P. Renard, Physica B+C 136B, 329 (1986); 
M. T. Hutchings, P. Day, E. Janke, and R. Pynn, J. Magn. 
Magn. Mater. 54-57, 673 (1986); S. T. Bramwell, M. T. 
Hutchings J. Norman, R. Pynn, and P. Day, J. de Phys. 
49, C8-1435 (1988); D. G. Wiesler, H. Zabel, and S. M. 
Shapiro, Physica B 156-7, 292 (1989); D. G. Wiesler, H. 
Zabel, and S. M. Shapiro, Z. Physik B 93, 277 (1994); 



L. P. Regnault, C. Lartigue, J. F. Legrand, B. Farago, J. 
Rossat-Mignod, and J. Y. Henry, Physica B 156-7, 298 
(1989). 

17. P. Gaveau, J. P. Boucher, L. P. Regnault, and Y. Henry, 
J. Appl. Phys. 69, 6228 (1991). 

18. D.L.Huber, Phys. Rev. B26, 3758 (1982). 

19. The sign of our damping term differs from Q|l^,^] be- 
cause we work with spins while these authors deal with 
magnetisations . 

20. S.Iida, J. Phys. Chem. Solids 24, 625 (1963). 

21. T. Kamppeter, F. G. Mertens, A. Sanchez, N. Gr0nbech- 
Jensen, A. R. Bishop, and F. Dommguez-Adame, in "The- 
ory of Spin Lattices and Lattice Gauge Models," eds. M. 
L. Ristig, K. A. Gernoth and J. W. Clark. Springer- Verlag, 
Berlin (1997). 

22. The correlation between different components of F st can 
also be calculated and turns out to be zero. 

23. S.Miyashita, H. Nishimori, A.Kuroda and M.Suzuki, 
Progress of Theor. Phys. 60, 1669 (1978). 

24. H.-J. Schnitzer, Zur Dynamik kollektiver Anregungen in 
Hamiltonschen Systemen, Ph.D.-thesis, University of 
Bayreuth (1996). 

25. N. Gr0nbech- Jensen and S. Doniach, J. Comp. Chem. 15, 
997 (1994). 

26. F. G. Mertens, G. M. Wysin. A. R. Volkel and A. R. Bishop, 

27. Strictly speaking, Eq. ( |48| ) must be solved in polar coor- 
dinates which leads to shifts of Wi,2 by ±Wo, respectively 
Jllj . The different signs appear due to the phase differences 
±7r/2, below Eq. (38). 

28. Yu. Gaididei, T. Kamppeter, F. G. Mertens and A. R. 
Bishop, submitted (August 1998). 

29. Yu. Gaididei, T. Kamppeter, F. G. Mertens and A. R. 
Bishop, in preparation. 

30. D. A. Dimitrov and G. M. Wysin, Phys. Rev. B 53, 8539 
(1996). 

31. J. E. R. Costa, B. V. Costa, and D. P. Landau, Phys. Rev. 
B 57, 11510 (1998). 

32. D. A. Dimitrov and G. M. Wysin, preprint (March 1998). 

33. J. Schliemann, F. G. Mertens, and H. Frahn, in prepara- 
tion. 



